home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / matrix.lha / Matrix / src / cosft.C next >
C/C++ Source or Header  |  1993-02-25  |  1KB  |  67 lines
  1. #include <math.h>
  2. #ifndef PI
  3. #define PI M_PI
  4. #endif
  5.  
  6. void cosft(double y[], int n, int isign)
  7.  
  8. /*
  9.   Calculates the cosine transform of a set y[1..n] of real-valued data points.
  10.   The transformed data replaces the original data in array y.
  11.   n must be a power of 2.
  12.   Set isign to +1 for a transform and to -1 for an inverse transform.
  13.   But the inverse transform array should be multiplied by 2/n.
  14. */
  15.  
  16. {
  17.   int     j, m, n2;
  18.   double enf0, even, odd, sum, sume, sumo, y1, y2;
  19.   double theta, wi = 0.0, wr = 1.0, wpi, wpr, wtemp;
  20.   void     realft(double*, int, int);
  21.  
  22.   theta    =  PI/(double) n;
  23.   wtemp    =  sin(0.5*theta);
  24.   wpr    = -2.0*wtemp*wtemp;
  25.   wpi    =  sin(theta);
  26.   sum    =  y[1];
  27.   m    =  n >> 1;
  28.   n2    =  n + 2;
  29.   for (j = 2; j <= m; j++)
  30.     {
  31.       wr     = (wtemp = wr)*wpr - wi*wpi + wr;
  32.       wi     =  wi*wpr + wtemp*wpi + wi;
  33.       y1     =  0.5*(y[j] + y[n2-j]);
  34.       y2     = (y[j] - y[n2-j]);
  35.       y[j]     =  y1 - wi*y2;
  36.       y[n2-j]     =  y1 + wi*y2;
  37.       sum    +=  wr*y2;
  38.     };
  39.         realft(y, m, 1);
  40.   y[2]     =  sum;
  41.   for (j = 4; j <= n; j += 2)
  42.     {
  43.       sum    +=  y[j];
  44.       y[j]     =  sum;
  45.     };
  46.   if (isign  ==  -1)
  47.     {
  48.       even     =  y[1];
  49.       odd     =  y[2];
  50.       for (j = 3; j <= n-1; j += 2)
  51.     {
  52.       even    +=  y[j];
  53.       odd    +=  y[j+1];
  54.     };
  55.       enf0     =  2.0*(even - odd);
  56.       sumo     =  y[1] - enf0;
  57.       sume     = (2.0*odd/n) - sumo;
  58.       y[1]     =  0.5*enf0;
  59.       y[2]    -=  sume;
  60.       for (j = 3; j <= n-1; j += 2)
  61.     {
  62.       y[j]        -=  sumo;
  63.       y[j+1]    -=  sume;
  64.     };
  65.     };
  66. }
  67.